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^ ■ Abstract 

D ■ 

A numerical simulation of fluid flows in a Laval nozzle is performed to observe the formation of 
an acoustic black hole and the classical counterpart to Hawking radiation under a realistic setting 
of the laboratory experiment. We aim to construct a practical procedure of the data analysis 
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\q ■ to extract the classical counterpart to Hawking radiation from experimental data. Following our 

so : 

procedure, we determine the surface gravity of the acoustic black hole from obtained numerical data. 
Some noteworthy points in analyzing the experimental data are clarified through our numerical 
■ simulation. 
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I. INTRODUCTION 

The existence of the black hole horizon is essential to the Hawking radiation ]I Q> 0, ^ • 
The horizon causes an extremely large redshift on an outgoing wave of a matter field during 
propagating from a vicinity of the horizon to the asymptotically flat region. This redshift 
results in the Planckian distribution of quantum mechanically created particles at the future 
null infinity. Although the Hawking radiation has not yet been confirmed observationally, 
it has been pointed out that the similar phenomenon appears on phonons in a transonic 
flow of a fluid. This theoretical phenomenon in the fluid is called the sonic analogue of 
Hawking radiation, and the transonic fluid flow on which the sonic analogue of Hawking 
rad.atiou occurs is caHcd the acoustic b ,ae k hole flflflflflUHQ. lu the trausouic 
flow, the boundary between a subsonic and a supersonic region is the sonic point. The sonic 
point corresponds to the black hole horizon and called the sonic horizon. The region apart 
form the sonic point in the subsonic region corresponds to the asymptotically flat region 
in a black hole spacetime. The sound wave which propagates from a vicinity of the sonic 
point to the subsonic region receives the extremely large redshift and the sonic analogue of 
Hawking radiation appears on phonons. The Hawking temperature Tjj of the acoustic black 
hole is given by the gradient of the fluid velocity at the sonic point x = x s Q]: 
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where c s and v are the sound velocity and the fluid velocity, respectively. For an ordinary 
system with typical size ~ 1 m and c s ~ 340 m/s, we obtain T# ~ 4 x 10~ 10 K. Thus, in 
a practical experiment in a laboratory, it is very difficult to detect the sonic analogue of 
Hawking radiation with such an extremely low temperature. 

While it is very hard to detect the sonic analogue of Hawking radiation on phonons, 
which is a quantum origin, a classical sound wave with a sufficiently large amplitude can 
be available for detecting the effect of the acoustic black hole on a propagation of sound 
waves in practical experiments. Indeed, the thermal nature of Hawking radiation is not 
a quantum origin; it comes from a large redshift on the wave propagating from a vicinity 
of the horizon to the asymptotically flat region, and has nothing to do with the quantum 
effect (see Appendix). Therefore, if we can construct a classical quantity from sound waves 
which shows the thermal distribution, then such a quantity can be interpreted as a classical 
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counterpart to Hawking radiation. A candidate for the classical counterpart has already 



14 1 by using a Fourier component of an outgoing wave 



been introduced in the references |13l 
from a vicinity of the horizon. 

By the way, since a realistic fluid is composed of molecules, the effect of the molecule size 
discrepancy of the fluid becomes one of the important issues of the acoustic black hole. This 
effect appears in the modified dispersion relation at high energies of phonons. However, in 
the expansion by the molecule size discrepancy, the zeroth order form of the quantum ex- 
pectation value of the number operator (0| N |0) of phonons results in the thermal spectrum 
at least not so high energy region. That is, the effect of the molecule size discrepancy is 
not dominant in the sonic analogue of Hawking radiation. Therefore, before proceeding to 
the issue of the molecule size discrepancy, it is necessary to detect the ordinary Hawking 
radiation (thermal spectrum) by a real experiment in laboratory, and the molecule size dis- 
crepancy is not in the scope of this paper. Hence we concentrate on the construction of the 
practical procedure to extract the evidence of the Hawking radiation (thermal spectrum) 
from the experimental data. To do so, we make use of the classical counterpart to Hawking 
radiation. 

We should emphasize here that the Hawking temperature (JIJ includes the Planck constant 
h and we can not consider the temperature Th in discussing the "classical" counterpart 
to Hawking radiation. However we can consider the surface gravity k of the black hole 
horizon which is a purely classical quantity and relates to the Hawking temperature as 
ksTn = hn/2n. Hence, for the acoustic black hole, we consider the surface gravity Kh of 
the sonic horizon instead of the Hawking temperature, 

d fv ^ (2) 



K H = C s — . 
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In this paper, we aim to construct a practical procedure of the data analysis to obtain 
the classical counterpart to Hawking radiation in a realistic setting of the acoustic black 
hole. To do so, we perform a numerical simulation of a fluid flow in a practical experimental 
setting, and demonstrate that our procedure works well to detect the classical counterpart 
to Hawking radiation. 

We organize the paper as follows. In the section |Hj we review the acoustic black hole in 
a Laval nozzle and define the classical counterpart to Hawking radiation. The section |^ 
is devoted to the introduction of our numerical method to simulate a transonic flow, and 



to the construction of the practical procedure of the data analysis. In the section IIV1 our 
numerical results are presented and we determine the surface gravity of the sonic horizon. 
Finally the summary and conclusion are in the section 

II. ACOUSTIC BLACK HOLE WITH A LAVAL NOZZLE 

We consider an acoustic black hole with a fluid in a Laval nozzle as proposed in Q| . The 
Laval nozzle, as shown in the right panel in FIG. has a throat where the cross section 
of the nozzle becomes minimum. The fluid is accelerated from the up stream to the down 
stream and the flow has the sonic point at the throat when a sufficiently large velocity is 
given at the inlet of the nozzle. Even if we prepare an initial flow which has no sonic point 
(e.g. the lower part of the right panel in FIG.Q), the flow can settle down to a stationary 
transonic flow (e.g. the upper part of the right panel in FIG.Q) with appropriate boundary 
conditions at the inlet and the outlet of the nozzle. This is the sonic analogue model which 
corresponds to the gravitational collapse. 

Here we note that, if a stationary transonic flow is prepared from the beginning (which 
corresponds to the eternal black hole), we can not expect to obtain the classical counterpart 
to the Hawking radiation; for the eternal case, the tunneling of phonons across the sonic 
point can result in the Hawking radiation and this is the purely quantum effect. Hence, 
when we are interested in the "classical" counterpart to Hawking radiation, it is necessary 
to consider the situation of the sonic point formation in course of the dynamical evolution of 
the fluid flow . As seen below, the sound wave which is prepared on the fluid flow before the 
formation of the sonic point becomes the classical counterpart to the quantum fluctuation 
which causes the classical counterpart to Hawking radiation. 
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FIG. 1: Propagation of waves in a spacetime of gravitational collapse and in a Laval nozzle. 
A. Basic equations 

We consider a perfect fluid and treat the flow in a Laval nozzle as one dimensional for 
simplicity. The basic equations are the mass conservation equation and the Euler equation: 



d t p+-d x (pvS) = 0, 
p (d t v + vd x v) = -d x p, 



(3a) 
(3b) 



where p(t,x) is the mass density, v(t, x) is the fluid velocity, p(t,x) is the pressure of the 
fluid and S(x) is the cross section of the Laval nozzle. We assume the adiabatic ideal gas 
type equation of state p oc p 7 where 7 is the adiabatic index. The sound velocity c s is given 
by 

(4) 
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For a stationary background flow, we can obtain the sound velocity c s and the cross section 
S of the nozzle as a function of the Mach number M = v/c s of the flow: 

-1/9 , - 7+1 

M- fM 2 + -Z T \ a «-v 

' ^ ' (5) 
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where quantities with subscript "in" represent the values at the inlet of the nozzle. The 
spatial derivative of the Mach number at the sonic point x = x s is 



dM 



dx 



± 
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fd 2 S\ 


4 S 


\dx 2 J 



(6) 



and this value relates to the surface gravity of the sonic horizon by Eq. (J2J) . Figure 121 shows 
the structure of the stationary flow for 7 = 7/5 and S(x) = 2 — cos(7rx). Here we set that 
the inlet is at x — 1, the throat at x = and the outlet at x — — 1. That is, the fluid flows 
downward from x = 1 to x = —1, then M < and FIG. |2 shows the absolute value of M. 
The boundary condition is given by the Mach number M in at the inlet of the nozzle, and 
each line in FIG. E] corresponds to a different value of M- m . The blue line which passes the 
sonic point (x,M) = (0, 1) is given by M- m = M* pa —0.197 and represents the stationary 
transonic flow. For the case q n = 1, we find 



dM 



dx 

and the surface gravity given by Eq. is 



±\/|t, c,U„ = (3IM.I) 1 / 8 , (7) 



^ = gg/„ 0.355. (8) 
2tt 2tt V 5 W 

Here it should be emphasized that, because the cross section of the nozzle is set S(x) = 
2 — cos(7rs) with — 1 < x < 1, the spatial scale is normalized by L/2, where L is the length 
of the nozzle, and that the temporal scale is normalized by L/(2c- m ), where we set c m = 1. 

We have not considered the perturbation of the fluid flow so far. In the following sections, 
the perturbation (sound wave) is introduced on the fluid flow, and the classical counterpart 
to Hawking radiation is defined. Then in the section IIV1 we compare this surface gravity 
(JHJ) with the surface gravity of the sonic horizon which is obtained through the observation 
of the classical counterpart to Hawking radiation in our numerical simulation. 
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FIG. 2: Absolute value of the Mach number as a function of a spatial coordinate x. Each line 
represents the stationary solution of the flow with different boundary value Mj„. The blue lines 
represent transonic flows. 

B. Classical counterpart to Hawking radiation 

As originally shown in the reference jsj], the perturbation of the velocity potential for 
the transonic fluid flow obeys the same equation as a massless free scalar field in a black 
hole spacetime. Therefore the classical counterpart to Hawking radiation of the acoustic 
black hole in the Laval nozzle should be observed by the sound wave on the transonic fluid 
flow. The most important sound wave is the one which starts to propagate against the 
stream from the outlet of the nozzle before the formation of the sonic point (the lower part 
of the right panel in FIG. and passes through the throat just before the moment of the 
sonic point formation to reach the inlet of the nozzle (the upper part of the right panel in 
FIG.Q). This sound wave receives the extremely large redshift which will cause the classical 
counterpart to Hawking radiation. 

When the fluid flow is irrotational and we introduce the velocity potential $ which relates 
with the fluid velocity as v — c\.$, the evolution equation of the perturbation of $ is 
obtained from Eqs. ©, 

[d t + d x v] [d t + vd x ] (j> - -Ld x (S pd x( j)) = 0. (9) 

Sp \cj J Sp 
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This corresponds to the Klein-Gordon equation □</> = with the acoustic metric 

ds 2 = (^^j [-( c l ~ v2 ) df2 ~ 2v dt dx + dx2 + d V 2 + dz<2 ] > ( 10 ) 

where the location of the sonic horizon (sonic point) x = x s is given by c s (x s ) 2 = v(x s ) 2 . 
Here we should note that, because we consider the one dimensional flow along x-axis as 
mentioned at the beginning of the previous section III Al the two dimensional section of 
(y, z) coordinate is omitted in obtaining Eq. Q from 0(f) = 0. In the following analysis, we 
omit the two dimensional section of (y, z) coordinate. 

In order to describe the extremely large redshift which causes the classical counterpart 
to Hawking radiation, we introduce the following null coordinates for the first, 

/dx f dx , . 

— — , w = t+ , 11 
c s + v J c s -v 

and the acoustic metric becomes ds 2 = (S p / c s ) (c 2 — v 2 ) du dw , which has a coordinate 

singularity at x = x s . In order to eliminate the coordinate singularity, we introduce the new 

null coordinates 

U = --^—exp[-(c' s0 + v' )u], W = -^—exp[(c' s0 + v' )w] J (12) 

where ' = d x and quantities with subscript denotes values at the sonic point x = x s . Then 
the acoustic metric and the evolution equation (jHJ) of the perturbation become 



,/.s- (^) 







f c s dx 


exp 


-2(c' s0 + ^ 











dUdW = -FdUdW, (13) 



□0=1^0 = 0. (14) 

The form of this metric near the sonic horizon is ds 2 ~ — 2SoPo(c' s o + v' ) dU dW, where 
c 2 — v 2 ~ 2c s0 (c' s0 + v' )x is used. This denotes explicitly that the coordinate (U, W) does 
not have coordinate singularities and corresponds to the Kruskal-Szekeres coordinate of the 
Schwarzschild spacetime. Here it should be recalled that, in the spacetime of the gravita- 
tional collapse (the left panel in FIG.HJ), the Kruskal-Szekeres coordinate is appropriate to 
describe the rest observer at the spacetime region before the formation of the black hole . 
Hence the outgoing normal mode of the sound wave which corresponds to the zero point 
fluctuation of the quantized matter field before the formation of the black hole (i.e. the 
mode function on the flat spacetime) is given by 

r\U) = A^e'^ u , (15) 
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where ujq is the initial frequency and As is the amplitude of this mode. In the original 
coordinate (t,x), this mode behaves as 

dx 



cf> m (U(u)) = cp m (t,x) = A 6 exp 



ZUJQ 



7 exp -{c' s0 + v' ) [t 



(16) 



. C> S + V ' I 1 ™ \ J Cs + V, 

This yields the temporal wave form at the observation point x = x D b s in the asymptotic 
region: 

iu 



out (t) =<p™(t,x ohB ) = A 6 exp 



Fourier component of this mode is defined by 



7 exp -(c' s0 + v' ) t 



^ obs 



dx 



c s + v 



tout 



and the power spectrum is 



4 2 



dtcj) out (t)e 



iujt 



(17) 



(18) 



where the quantity ujh is defined by 



Cso_d_f v 
2n dx \ c s 



— OO < LU < +OOj, 



(19) 



(20) 



Here it should be noted that, because the existence of the sonic point is assumed in the 
above discussion, the form of P{oS) of Eq. (|T9*|) must be obtained after the formation of the 
sonic horizon in the fluid flow. 

One may think it is strange that the initial frequency ujq does not appear in the power 
spectrum (119)1 . It is the fact that the initial frequency ujq does appear in the power spectrum 
in the context of the quantum field theory, because the amplitude is determined by the 
normalization condition for the mode functions to be A^ oc 1/y/uJo. However, in the calcula- 
tion for the classical counterpart to Hawking radiation, the initial frequency c^o appears only 
in the phase of the Fourier component 0° ut and the power P{uS) does not explicitly depend 
on Uq. 

The power spectrum (|19|) has the Planckian distribution for uo < and this is the classical 
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141 ] . For the quantum Hawking radiation, 



counterpart to the quantum Hawking radiation 
the magnitude of the power is determined by the zero point oscillation of the quantized field. 
However, for the classical counterpart to Hawking radiation, the magnitude of the power 
depends on the amplitude As of the input signal and it is expected that we can detect 



the thermal distribution of the power spectrum for the sound wave of a sufficiently large 
amplitude in practical experiments. 

As long as concerning one dimensional fluid flow, the amplitude of the sound wave does not 
decrease during propagating from the outlet to the inlet of the nozzle. The one dimensional 
transonic fluid flow corresponds to the two dimensional black hole spacetime in which no 
curvature scattering occurs. The theoretically expected power (JT5j) is equivalent to the 
formula of the Hawking radiation for black holes derived ignoring the curvature scattering. 

Our purpose is to observe this thermal spectrum via the numerical simulation of a tran- 
sonic flow in the Laval nozzle, and to construct the practical procedure of the data analysis. 
Because the observable sound wave is a real valued one, we prepare the following real valued 
input mode at the outlet of the nozzle before the formation of the sonic horizon: 

in ^(t) = A^ C os(u t + 9) = ~ [e- %e ^{U{t,x c )) + e w ^*(U(t,x c ))\ , (21) 

where 9 is the initial phase of the wave, <p m (U) is given by Eq. ([15)1 . and x = x e is the 
location of the outlet of the nozzle. Here we should note that, although the outlet of the 
nozzle is the "exit" of the background fluid flow, the sound wave can propagate against the 
background flow from the outlet to the inlet of the nozzle. Then, since the amplitude of 
the sound wave on one dimensional fluid flow does not decrease, the input mode of Eq. (j21j) 
causes the following sound wave in the fluid flow, 

^(t,x) = i [e- w r(U(t,x))+e w r*(U(t,x))] . (22) 

The effect of the redshift due to the formation of the acoustic black hole gradually appears on 
the observed sound wave. The temporal evolution of the observed sound wave is affected by 
the formation of the acoustic black hole. That is, the information of the classical counterpart 
to Hawking radiation is encoded in the observed sound wave. Therefore, using the output 
form (|17|). the resulting output signal after the formation of the sonic horizon is given by 

0°ut(0)^ = I [ e -i00out^ + gi^out*^] ; (23) 

and referring Eq. (|18|) . the Fourier component of this output signal is given by 

C W = \ (^Vr + zUT) ■ (24) 

Hence, because the initial phase 9 appears in 0^ ut ^' ) , we can not obtain the pure Fourier 
component 0° ut by observing once the real valued output signal (j) out ^ e \t). In order to 
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retrieve the pure Fourier component </>° ut from the mixed Fourier component <jh^ , we have 
to superpose two output modes out ( 6 ') and out ( 6l + 7r / 2 ) ) and obtain 

<PT = e ie (0° ut(6) + i<tC {6+7T/2) ) ■ (25) 

It is obvious that the power spectrum of this Fourier component gives the same form as 
Eq. (fTTJj) . Therefore, in order to obtain the Planckian distribution in the power spectrum 
P{w) of the observed sound wave from real valued input signals, we have to combine at least 
two output modes of which the input phases differ by n/2. 



III. NUMERICAL SIMULATION OF TRANSONIC FLOWS IN A LAVAL NOZ- 
ZLE 

Our numerical simulation is designed to generate a one dimensional transonic fluid flow 
from the initial configuration with no sonic point. We prepare the input signal at the outlet 
of the nozzle and observe the sound wave at the inlet of the nozzle. At the same time, the 
practical procedure of the data analysis is constructed. 



A. Basic equations for numerical simulations 

For the numerical calculation, we rewrite the fluid Eqs. Q using c s and v, 

-d t c s H ^—vd x c s + c s d x v + ( %^ J c s v = 0, (26a) 



7 — 1 y o 

2 

d t v + vd x v H -c s d x c s = 0. (26b) 

7-1 

This set of equations is not suitable for a numerical simulation of wave propagation in a 
transonic flow. We transform these equations to the advection form. We introduce the 
Riemann invariants J±(t,x) as follows 

J ± = v±^-. (27) 
7 - 1 



Then the basic equations become 

9 J. = Jl^ltt 

S 

%>• 

~s 



d t J + + V + d x J + = -1-1 (^f) (J* - J 2 ) , (28a) 



cV- - = +1-1 (^f) (J\ - J!) , (28b) 
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where 



V + = c s + v = ^ J+ - 1^ J_, K_ = c s - « = ^ J + - ^±1 J_. (29) 



The left hand side of Eqs. ()28|) are of the advection form. That is, if d x S = and d t V± = 0, 
J + propagates upward along t — J dx/V + = const, and J_ propagates downward along 
t + J dx/V- = const. This form of the equations is suitable to treat the propagation of 
waves numerically. 

In the sonic analogue of Hawking radiation, the perturbation of the velocity potential 
(sound wave) corresponds to the scalar field in a black hole spacetime as shown by Eq. Q. 
Therefore, in a practical experiment of an acoustic black hole, we should calculate the 
velocity potential by integrating the observed velocity of the fluid: 

r 2-obs 

^(t,x ohs ) = I dxv(t,x), (30) 

J xo 

where x b s is the observation point of sound waves and xq defines the origin of the velocity 
potential. However, since our experimental setting is designed to observe the sound wave at 
a fixed spatial point x = x ^ s , the observational data is a temporal sequence of v(t, x bs) and 
we must devise the other method to obtain the velocity potential at the observation point. 
In the upstream subsonic region of the transonic flow where the effects of the sonic horizon 
is negligible and the flow is stationary, the sound wave propagates along the "null" direction 
as indicated by Eq. (JTTJ) 

x 

t — ; = const. (31) 

and we can obtain the value of the velocity potential at x Q b s by integrating the velocity with 
respect to time 

®(t,x ohs ) = (c s +v)\ Xoh I dtv(t,x ohs ), (32) 

° 8 J t 

where t defines the origin of the potential. We use this formula to evaluate the velocity 
potential at the observation point. 

B. Experimental setting and numerical method 

We use the Laval nozzle with the cross section of the following form, 

S(x) = 2 - cos (ttx) (-1<x<1). (33) 
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We make the fluid flow against x-axis and assume v < 0. The inlet of the nozzle is at x = 1, 
the throat at x = and the outlet at x = — 1. The input wave is prepared at the outlet 
of the nozzle and the wave propagates against the flow from x = — 1 to x = 1. After the 
flow settles down to a stationary transonic flow, since d x S\ =1 = at the inlet, the set of 
Eqs. (I28|) gives d x v\ 1 = at the inlet. This means that the inlet of the nozzle corresponds 
to the asymptotically flat region in the gravitational collapse spacetime. Hence we make the 
observation of the sound wave at the inlet and set x D b s = 1. 

We prepare a flow with a homogeneous velocity distribution as the initial configuration, 

v(0, x) = Vi = const. < 0, c s (0, x) = 1. (34) 

This initial configuration of flow has no sonic point and corresponds to the flat spacetime 
region before the formation of a black hole. Since our evolution equations (J28|) are of the 
advection form, it is enough to set the boundary conditions for J + at the outlet and J_ at 
the inlet: 

J+(t,-l) = J+(0,-l) + Aj cos (uj t + 9) , J_(t, 1) = J_(0,1) (35) 

where Aj is a constant that represents the amplitude of the perturbation of J + at x — — 1, 
and the values J+(0, —1) and J_(0, 1) are determined consistent with Eq. (|2"7jl and Eq. (J33J). 
These boundary conditions mean that sound waves of constant amplitude continue to emerge 
from the outlet toward the inlet of the nozzle all the time of the numerical simulation. Here 
recall that the perturbation <fi of the velocity potential <£> plays the role of the scalar field in 
the ordinary Hawking radiation. The second term in J+(t, — 1) at the outlet Aj cos (uot + 6) 
causes perturbations of the velocity potential and results in the classical counterpart of 
Hawking radiation at the inlet of the nozzle. Because the amplitude of the sound wave in 
one dimensional fluid flow does not decrease, the amplitude of the input mode of the velocity 
potential (|2"T]) is given by 

where c s \ and vi are respectively the sound velocity and the fluid velocity at x = 1, and the 
relation of the perturbations <j>(t, x = —1) = 6&(t, —1) = (c s i + V\) J 1 dtSJ + (t, — 1)/2 at the 
outlet is used (see Eqs. f)27p and 

In order to extract the perturbation part from the "full" velocity potential $ of Eq. (|32|) 
which includes the "background" flow of the fluid, we should prepare three input modes 
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with different initial phases 9\ = —71/2,62 = and #3 = tt/2. If the amplitude Aj of 
the perturbation part of J+ is small enough, then the backreaction effects in the velocity 
potentials $ out ( 6 '*) [i = 1,2,3) at the observation point (x = 1) are negligible and we can 
subtract the common background contribution by taking the difference among them. Hence 
we obtain the Fourier components of the perturbation (sound wave) at the observation point 
as follows 

0out(12) _ (jjout(0i) _ (jj out (0 2 ) _ 0out(0i) _ 0out(0 2 )^ (37a) 
0out(23) _ (j>out(0 2 ) _ ^out(0 3 ) _ 0out(0 2 ) _ out ( e 3^ (37b) 

where $2; Ut is the Fourier component of $ out ( *). These Fourier components correspond 
to the real valued output signal of Eq. (|24|) . Here we note that, according to the boundary 
condition (|3*5)l . the perturbation parts of J + which produce the output signals of Eqs. (jHZj) 
are given by 



5J^ 2) = V2Aj cos (u t + ^pj , 5 j| 23) = V2Aj cos (u t + - 



(38) 



where we used cos(u;ot) — cos(u;ot ± 7r/2) = ^/2cos(ujot ± 7r/4). Hence we obtain the pure 

Fourier component 0° ut of the output sound wave by Eq. (|23j) . and its power spectrum is 
given by 

| out|2 = l| C ut ( 23) + . out(12)|^ (39) 

where the factor 1/2 is introduced to eliminate the factor \/2 which appears in the amplitude 
\/2Aj of the input signals (|3*5j) . When we carry out a realistic experiment of the acoustic 
black hole, the power spectrum gives the classical counterpart to Hawking radiation 
and it has to be compared with the theoretical form (|19j) . Hence, it is recognized that we 
have to perform at least three independent experiments with three different input phases 
and observe three independent output signals to obtain the classical counterpart to Hawking 
radiation. 

With our setting of the fluid flow in the Laval nozzle, we carry out the numerical simu- 
lation by the finite difference method. We use the Cubic-Interpolated Pseudoparticle (CIP) 
method [15] for the interpolation between neighboring spatial mesh points. This method 
enables us to treat the shock which appears after the formation of the sonic point. We note 
here that the shock arises in the supersonic region and it does not affect the subsonic region 
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where we observe the output signal. The procedure of the numerical simulation to observe 
the thermal power spectrum (fT^j) of the perturbation of the velocity potential is as follows: 

Step 1: Generate the transonic fluid flow three times with three different initial phase 
values, Q\ = — n/2, 62 = and #3 = tt/2. Then, using Eq. (|32|). obtain the full velocity 
potentials $ out (^) (i = 1,2,3). 

Step 2: Calculate the quantities out ( 12 ) = $ out ( 9 i) anc j 0out(23) = $out(e 2 ) _ $ out(0 3 )_ 
Then compute the Fourier components of them, 0^ ut( - 12 ' ) and 0£, ut< ' 23 ' ) . 

Step 3: Calculate the power spectrum 

PoM = l\rj t{23) + *<f>T il2) \ 2 - (40) 

This is the observationally obtained power spectrum for the perturbation of the ve- 
locity potential. 

Step 4: Then compare this power spectrum with the theoretically expected one 
p ,s f (csi + vi) . \ 2 1 e»/»" __k h c s0 fv\ 

^theory (w) = — Aj -— , U H = — = — - , 41 

V 2uq J oJu^e, 1 H — 1 27r 2n \c s ' 



x — ,r * 



and determine the surface gravity kh of the sonic horizon. Here note that the quantity 
Cj in the amplitude of P t heor y is the effective input frequency given by the following 
discussion. 

We must note here that there arises an extra redshift effect on the sound wave. In our 
numerical simulation, the initial fluid flow has homogeneous velocity distribution. Then the 
system starts to evolve dynamically and settles down to the stationary transonic fluid flow 
finally. Therefore, the observed frequency shifts due to the evolution of the "background" 
flow. According to the null coordinate (fTTj) . the outgoing wave with frequency Uq is 

dx 



(pit, x — 1) oc exp ^— iojq (t — j 



c s + v 

and the effective frequency of the wave at x = 1 is given by 



(42) 



The second term in Eq. (J43|) gives a negative contribution in our non-stationary setting of 
numerical experiment and the observed frequency uq at x = 1 becomes smaller than ujq. 
The theoretically expected form of the power spectrum in our setting should be given by 
replacing Uo to lj in the formula (fTTJj) . 
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IV. NUMERICAL RESULTS 

We use the parameters 7 = 7/5, Aj = 1.0 x 1CT 6 and ojq = 100. The number of spatial 
mesh points are 10001 and the size of one mesh becomes Ax = 2/10001. The size of one 
temporal step is set At = Ax/10, and we calculate 900000 steps in time. The numerical 
simulation runs in the temporal range t = ~ 18. The spatial scale is normalized by L/2, 
where L is the length of the nozzle. The temporal scale is normalized by L/2c s (t = 0), where 
we set c s (t = 0) = 1 as denoted by the initial condition (|51jl. 

A. no black hole case 



For the initial fluid velocity V{ = —0.19, we do not observe the formation of an acoustic 
black hole. Figure H3 shows the evolution of the Mach number. After t ~ 10, the Mach 
number at x = reaches a constant value and a stationary flow is realized. 




FIG. 3: The left panel shows the temporal evolution of the spatial distribution of the Mach number. 
The right panel shows the Mach number at x = 0. 
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FIG. 4: Spatial distribution of the Mach number at t = 18. 



The spatial distribution of the fluid velocity at t — 18 (FIG.EJ) coincides with the stationary 
solution (jSJ) with M in = —0.19 and no sonic point appears. Figure El shows the time derivative 
of the observed velocity at x — 1 which should disappear when the fluid flow becomes 
stationary. Until t ~ 4, the initial burst mode appears, which is peculiar to our initial and 
boundary conditions. Then, the output signal decays exponentially. For Aj ^ (with input 
perturbation), the output signal oscillates with a constant amplitude after the "background" 
flow settles down to the stationary flow. 




t t 

FIG. 5: The temporal evolution of v at x = 1. The left panel is Aj = case and the right panel 
is Aj / case. 

Figures El and [7| show the spacetime distribution of the velocity perturbation for A j ^ 
with ujq = 20. 1 In FIG.I3 to subtract the "background" flow, we have defined Sv = v^—v^ 

1 For FIG. El 171131 and El we use d the smaller value of the perturbation frequencty uj — 20 to visualize 
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where i>M denotes the fluid velocity with the input phase 0j. We can see that the sound wave 
propagates from x — — 1 to x — 1. Around the throat x = 0, the fluid velocity has larger 
value compared to the other region and the phase velocity of the outgoing wave becomes 
smaller. Thus the slope of the constant phase line becomes larger around the throat. 




-1 □ 1 



FIG. 6: Spacetime distribution of v for Vi = —0.19 with Aj ^ 0, ujq = 20. Colors are assigned 
according to the value of log \v\. 

the propagation of the wave in the spacetime diagram. 
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FIG. 7: Spacetime distribution of the velocity perturbation Sv = — v^ 1 ' for Vi = —0.19 with 
Aj 7^ 0,uo = 20. Colors are assigned according to the value of log \5v\. 



Figure |H1 shows the perturbation of the velocity potential at x = 1. We set t = 3 
in Eq. to obtain the velocity potential. The observed perturbation of the velocity 
potential oscillates with a constant amplitude. 
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FIG. 8: The temporal evolution of the velocity potential out ( 12 ) and out ( 23 ). 
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FIG. 9: The temporal evolution of the power spectrum P b s (u;) for the velocity potential. 



The power spectrum P Q b s (w) obtained from the perturbations are shown in FIG. El The 
frequency of the observed signal evolves in time. This is due to the non-stationarity of the 
background flow as denoted by Eq. Until the flow settles down to the stationary flow 

which is determined by the given boundary condition, the background flow evolves in time 
and causes the shift of the frequency of the observed perturbation. After a sufficiently long 
time has passed and the fluid flow becomes stationary, the observed frequency coincides with 
the input frequency ojq = 100 as expected by Eq. (J33J) . 



B. 



black hole formation case 



For the initial fluid velocity Vi = —0.2, we observe the formation of the acoustic black 
hole. Figure ITU1 shows the evolution of the Mach number. 
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FIG. 10: The left panel shows the temporal evolution of the spatial distribution of the Mach 
number. The right panel shows the Mach number at x = 0. 
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FIG. 11: Spatial distribution of the Mach number at t = 18. 

At t ~ 10, the Mach number at the throat x = reaches unity and the sonic point appears 
(formation of the acoustic black hole). At the same time, the discontinuity of the derivative 
of the Mach number appears in the supersonic region x < as shown in FIG. ^2 This 
discontinuity is due to the shock formed in the transonic flow. The formation of the shock 
is peculiar to our boundary condition which determines the Riemann invariants J± at the 
inlet and the outlet of the nozzle. However, since the shock occurs in the supersonic region 
after the formation of the sonic point, the effect of the shock never propagate into the 
subsonic region after the formation of the sonic point. Hence the shock never affect the 
sonic analogue of Hawking radiation, and we do not have to pay attention to the effect of 
the shock formation. 
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Figure IT2l is the time derivative of the fluid velocity at x = 1. After the burst mode t <~ 4, 
we can observe the quasi-normal oscillation of the acoustic black hole with a period ~ 4. 




FIG. 12: The temporal evolution of v at x = 1. The left panel is = case and the right panel 
is Aj y£ case. As the amplitude of the input perturbation is far smaller than the change of the 
background, we can not see the oscillation of the perturbation in the figure (right panel). 

The value of this period can be estimated as follows. In our simulation, the value of J_ 
at the inlet of the nozzle is fixed to be constant by the boundary condition. Hence after 
the formation of the sonic horizon, the boundary condition for the ingoing perturbation 8 J_ 
becomes 

5J_ = free : at the horizon, 5J_ = : at the inlet of the nozzle (44) 

For the outgoing perturbation 5J + , 

5J + = : at the horizon, 5J + = free : at the inlet of the nozzle (45) 

By these boundary condition, a quarter of the wavelength of the quasi-normal mode is equal 
to a half length of the Laval nozzle and the wavelength of the quasi-normal mode becomes 
4 in the present case. Thus, we obtain the period 4 of the quasi-normal mode by assuming 
c s = 1. 

The spacetime distribution of the time derivative of the fluid velocity and the perturbation 
is shown in FIG. UHland FIG. El respectively. 



22 



FIG. 13: Spacetime distribution of v for Vi = —0.2 with Aj ^ 0,wo = 20. Colors are assigned 
according to the value of log \ v\. 
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FIG. 14: Spacetime distribution of the velocity perturbation 5v = v&> — v^' for Vi = —0.2 with 
Aj ^ 0,uo = 20. Colors are assigned according to the value of log \5v\. 

The time evolution of the observed perturbation of the velocity potential is shown in 
FIG. ED After the formation of the sonic horizon t ~ 10, the frequency of the observed 
perturbation gradually decreases to zero. 




2.5 5 7.5 10 12.5 15 17.5 2.5 5 7.5 10 12.5 15 17.5 



t t 

FIG. 15: The temporal evolution of the velocity potential out ( 12 ) and (/> out ( 32 ). 

The time evolution of the power spectrum of the observed signal is shown in FIG. El 
FIG. [T7|and FIG. [TSJ The effective frequency Cjq of the input perturbation observed at the 
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inlet of the nozzle is 74.6. Figure fTBI shows that, after t ~ 10, the power spectrum of the 
observed perturbation spreads out toward the low uj range by the redshift effect due to the 
formation of the sonic horizon, and the divergence of the power P ^ s (uj) appears at uj = 0. 
This indicates that the observed spectrum approaches the theoretically expected form (J41j) 
which diverges at uj = as oc 1/uj. To remove this divergence, we plot uj 2 P ^, s (uj) in FIG. IT71 
and FIG. d 
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FIG. 16: The temporal evolution of the power spectrum P b s (w) for the output signal. 



25 



1.75x10" 
1.5x10" 
1.25x10" 



1x10 
7.5x10' 

5x10' 
2.5x10' 



-12 
■13 



■13 
■13 



t = 6 



1.75x10" 
1.5x10" 
1.25x10" 
1x10" 
7.5x10" 
5x10" 
2.5x10 



-75 -50 -25 

CO 



25 50 75 100 



-13 



t = 12 








J 


I 



-75 -50 -25 25 50 75 100 

CO 



1.75x10 
1.5x10 
1.25x10 
1x10 
7.5x10 
5x10 
2.5x10 



-12 
-12 



-12 
-12 



-13 
-13 



-13 



t = 14 



-75 -50 -25 

CO 



1.75x10" 
1.5x10" 
1.25x10' 
1x10" 
7.5x10' 
5x10" 
2.5x10' 



25 50 75 100 




-75 -50 -25 

CO 



25 50 75 100 



FIG. 17: The temporal evolution of the power spectrum u> 2 P \y S (uj) for the output signal. 
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FIG. 18: The temporal evolution of the power spectrum uj 2 P ^, s (lu) for the output signal in the 
range — 15 < u < 15. 
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C. Surface gravity of the sonic horizon 



Figure Uni shows the observed power spectrum a; 2 P bs(^) at £ = 18 given by Eq. (flTjj) and 

theory ( 



the theoretically expected power co> 2 P thcory (u;) given by Eq. (J4*T]) 




CO 




FIG. 19: Power spectrums uj 2 P ^, s {ui) in different ranges of ui. The solid lines are the observed 
power spectrums and the red dotted lines are theoretically expected power (|41|) with ujh = 0.360. 



We fits the theoretical power spectrum to the numerical one in the range — 10 < u < 10 
with the parameters c s \ = 1.00, V\ = — 0.20, Aj = 1.00 x 10~ 6 and Oq = 74.6, where the 
value of ujq is read off from the location of the highest peak of the power in FIG. El Then 
we obtain the surface gravity of the sonic horizon 

u H = — = 0.360 ±0.013. (46) 
2n 

This value is consistent with the theoretically expected value 0.355 given by Eq. (JHJ). 

Here we discuss about two points. The first is about the oscillation of the numerical power 
spectrum. This oscillation is due to the finite size Fourier transformation. This oscillation 
will disappear if the numerical calculation is carried out with larger value of lvq. 

The second point is the deviation of the numerical power from the theoretical one in the 
high u range. This deviation is due to the input sound wave given in the boundary condition 
(jHEJ). The input sound wave emerges from the outlet of the nozzle toward the inlet of the 
nozzle, and the observation is done at the inlet. Before the formation of the sonic point, 
the input sound wave comes from the outlet to the inlet without receiving a large redshift. 
Then, as the fluid flow evolves in time, the redshift on the observed sound wave becomes 
larger, and the observed frequency continues to decrease until the sonic point appears, as 
denoted by Eq. (|4*3*|) . Hence, although the input sound wave is a monochromatic wave at the 
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outlet of the nozzle (the starting point of the propagation), the observed sound wave at the 
inlet of the nozzle has a broad power spectrum after the formation of the sonic point because 
of the stimulated effect of the redshift before the formation of the sonic point. The peak of 
this broad spectrum appears obviously in FIG. El around u ~ 75. This broadness in the 
power spectrum is the origin of the deviation of the numerical power from the theoretical 
one in FIG. [03 The perfect thermal spectrum should be realized at the infinite future, 
and that the theoretical power spectrum ()41|) is the form evaluated at the infinite future. 
However we have to make the observation in a finite temporal interval. Therefore, when we 
plot the numerical result in an appropriately large range of a;, it is impossible to avoid some 
deviation of the numerical power from the theoretical one. If we carry out the numerical 
calculation longer and longer, a better agreement have to be obtained in FIG. [TTfl 

From above discussions, we conclude that FIG. IT^l shows the good agreement between the 
numerically obtained power and the theoretically expected form of the classical counterpart 
to Hawking radiation. 

V. SUMMARY AND CONCLUSION 

For the acoustic black hole, the classical counterpart to Hawking radiation is given by 
the power spectrum of the perturbation of the velocity potential of the fluid and detectable 
in practical experiments. To demonstrate its detectability, we performed the numerical sim- 
ulation of the acoustic black hole in the Laval nozzle and observed the classical counterpart 
to Hawking radiation. We obtained the good agreement of the numerically observed power 
spectrum with the theoretically expected one. 

Through our numerical simulation, we have obtained two noteworthy points for data 
analysis of the experiments of the acoustic black hole: the first one is that a single input 
wave can not give us necessary information of the classical counterpart to Hawking radiation. 
We need to carry out several independent observations of the sound waves with different 
initial phases to retrieve the thermal distribution of the power spectrum. The second one is 
that we can evaluate the velocity potential at the observation point by integrating the fluid 
velocity at that point with respect to time (see Eq. (|32|l). Therefore, it is not necessary to 
observe the sound wave at every spatial points in the fluid. 

In our calculation presented here, we have considered the sound wave with small ampli- 
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tude. However our numerical code is applicable beyond perturbation, in which the non-linear 
effect of the sound wave becomes important. By analyzing such a situation, we may be able 
to discuss the backreaction effect on the classical counterpart to Hawking radiation. As the 
next step, we are planning to extend our simulation of transonic flows in a Laval nozzle to 
include quantum effects. Then, we expect to obtain implications for quantum aspects of 
Hawking radiation using the sonic analogue model of black hole. 
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APPENDIX: CLASSICAL AND QUANTUM EFFECTS IN HAWKING RADIA- 
TION 

We briefly review the Hawking radiation in a spacetime of a gravitational collapse forming 
a Schwarzschild black hole, and clarify the distinction between the quantum effects and the 
classical effects in the occurrence of the Hawking radiation. For simplicity, we consider a 
massless free scalar field $ as a representative of the matter field, and set c = — 1. 

For the first, we start with the classical effects. Because the spacetime is dynamical, the 
positive frequency modes and the negative frequency modes of the scalar field are mixed as 
the system evolves in time. This mixing is represented by the Bogoliubov transformation 
between the positive frequency mode W of $ at the past null infinity and that mode 0^ at 
the future null infinity 



where 0* = is the negative frequency mode at the future null infinity. The Bogoliubov 
coefficients A and B are obtained by solving the classical wave equation □$ = with 
appropriate boundary conditions at the past and future null infinities and with the following 
relation, 



which comes from the normalization condition , (j)^) = 8 W1U2 with respect to the Klein- 
Gordon inner product of $ [2 J . This means that the derivation of A and B is purely classical. 

Next, in the classical framework, we consider the wave mode which propagates from the 
past null infinity to the future null infinity via a vicinity of the black hole horizon (see the left 
panel in FIG.Q). This wave mode is ingoing at the past null infinity and becomes outgoing 
at the future null infinity. The ingoing positive frequency mode 0^ at the past null infinity 
is given by 

(puijw) = exp (— iuw) , u> > , (A. 3) 

where w is the ingoing null coordinate appropriate for the rest observer at the past null 
infinity. When this wave mode <f> u (w) propagates to the future null infinity via a vicinity of 




(A.l) 




(A.2) 
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the horizon, it evolves to be (p^iwiu)) at the future null infinity under the geometrical optics 
approximation, where u is the outgoing null coordinate appropriate for the rest observer 
at the future null infinity. The function w(u) expresses the extremely large redshift which 
the mode 4>uj(w) receives during propagating from the past null infinity to the future null 
infinity. This redshift effect can be decomposed into two parts: one is the redshift during 
the propagation from the past null infinity to the vicinity of the horizon, and the other part 
is the redshift after passing through the vicinity of the horizon. The first contribution is not 
so large and the mixing of positive and negative frequency modes does not occur. However, 
the second contribution is large enough to make the Bogoliubov coefficient B w q non zero. 
By matching the null coordinates w and u along a null geodesic which connects the past and 
future null infinities via the vicinity of the horizon, the function w(u) is obtained [3,0: 

w(u) = —aexp (—ku) + 5, (A. 4) 

where k is the surface gravity of the black hole horizon, 5 is an arbitrary constant denoting 
the freedom of choosing the origin of w, and the constant a is determined by the first part of 
the redshift. The exponential form exp(— ku) in the Eq. (fOJ) comes from the second part of 
the redshift and implies that the wave length of the outgoing wave is exponentially stretched 
during propagating from a vicinity of the horizon to the future null infinity. That is, the 
wave <j)u(w(u)) is no longer a pure positive frequency mode but becomes a superposition of 
positive and negative frequency modes at the future null infinity. In order to calculate the 
superposition, we need the outgoing positive frequency mode W at the future null infinity: 

4>lu( u ) = —j= exp (— iuu) , oo > 0. (A. 5) 

Then, using the definition of the Bogoliubov transformation (jA.l|) . the wave 4> w (w(u)) at 
the future null infinity is decomposed as 

POO 

<f> u {w(u)) = / dw x A*;i0« 1 (w) + -Bwwifc 1 («) • ( A -6) 







The Bogoliubov coefficients are obtained using the inner product ((pwn&un) — A UlW2 and 



^,C 2 ) = --Bwjwa, and we find 



14 I 2 _ pW 2 /w fl I d p id | 2 , _ ( A 7\ 
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The square of the Bogoliubov coefficients do not depend on the constants a and 5, and the 
following relation holds: 

1-^0)1,— 0^2! l-^OJlO^I ' (■^■■^) 

All the above phenomena are the classical effects. 

Finally we proceed to the quantum effects. When <3> is quantized, the harmonic operators 
a w and with respect to the past mode are related to those fi^ and fit, with respect to 
the future mode W as follows P, y|: 

poo poo 

a u = du x (A UlU) a Ul + B* lW a* J , a* = / duJ i {K^ a U + B ^ a^) , (A.9) 
Jo Jo 

[ a m, a t 2 ] = &<W«u [a wi ,alj = ^<y wiwa . (A. 10) 

Then the number of particles at the future null infinity is obtained 

POO 

N u = (0\ala u \0) = h du^B^ 2 , (A.ll) 



JO 

where |0) is the vacuum state at the past null infinity, |0) = 0. Therefore, for the mode 
W which passes a vicinity of the horizon and propagates to the future null infinity, using 
Eqs. flA.7|) and an appropriate regularization method we obtain 

n u = (o\ at a. \o) = e ^ H _ 1 . (A.12) 

It is concluded that a black hole emits a thermal radiation of $ with the Hawking temper- 
ature T H = h un- 
it should be emphasized that, although the creation of particles is just the quantum effect, 
however the Planckian distribution (|A.12|) of the emitted particles is purely the classical 
effect due to Eqs. (jA.7|) . The thermal nature of the spectrum comes from the Bogoliubov 
coefficient |-B^ lW2 | 2 which has the Planckian distribution with respect to u%. That is, the 
thermal nature of the Hawking radiation comes from the extremely large redshift for the 
wave mode 0^ which passes the vicinity of the horizon. If we take the classical limit h — > 0, 
the number of created particles and the Hawking temperature become zero but \B\ 2 is not 
affected. 
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